model_run_summary

library(dplyr)
library(tidyr)
library(raster)
library(Matrix)
library(sparseMatrixStats)
library(Hmisc)
library(seegSDM)

library(ggplot2)

run_unique_name <- "MT_data_2"

# save this file before run. and ensure variable above is set correctly (before and after!)
rmarkdown::render('code_ruarai/R_interactive/summarise_model_run.Rmd',
                  output_file = paste0('reports/report_', run_unique_name, '.html'))

pred_files <- list.files(paste0("output/update/predictions/", run_unique_name, "/"), full.names=TRUE)

preds <- lapply(pred_files, function(x){
  pred_values <- readRDS(x)
  
  pred_values[is.na(pred_values)] <- 0
  
  return(Matrix::t(as(pred_values, "sparseMatrix")))
})
saveRDS(do.call(rbind, preds), 
        file = paste0("output/update/", run_unique_name , "_pred_matrix.Rds"),
        compress=FALSE)

full_pred_matrix <- readRDS(file = paste0("output/update/", run_unique_name, "_pred_matrix.Rds"))
blank_seasia <- raster('data/clean/raster/SEAsia_extent.grd')

cols <- colorRampPalette(c("#55843b", "#a4be79","#ffffbf", "#921d67", "#6c0043"))(100)

pred_means <- colMeans2(full_pred_matrix)
pred_sd <- colSds(full_pred_matrix)

pred_means[pred_means == 0] <- NA
pred_sd[pred_sd == 0] <- NA

plot(setValues(blank_seasia, pred_means), col = cols, zlim=c(0,1))
title("Predicted risk means", run_unique_name)

plot(setValues(blank_seasia, pred_sd))
title("Predicted risk standard deviation", run_unique_name)

risk_freya <- raster('data/clean/raster/SEAsia.tif')
diff_means <- pred_means - risk_freya

plot(diff_means, col=cols, zlim=c(-0.5,0.5))
title("Change in risk (vs FS)", paste0(run_unique_name, " - FS Risk"))

hist(diff_means, xlab='Change in risk from original FS risk map', main='')
title(paste0(run_unique_name, " - FS Risk"), run_unique_name)

# Change from arbitrary other dataset
comp_run <- "MT_data_1"

full_pred_matrix_comp <- readRDS(file = paste0("output/update/", comp_run, "_pred_matrix.Rds"))

comp_means <- colMeans2(full_pred_matrix_comp)

comp_diff_means <- pred_means - comp_means

plot(setValues(blank_seasia, comp_diff_means), col=cols, zlim=c(-0.5,0.5))
title(paste0("Change in risk (vs ", comp_run, ")"), paste0(run_unique_name, " - ", comp_run))

hist(comp_diff_means, xlab=paste0('Change in risk from ', run_unique_name), main='')
title(paste0(run_unique_name, " - ", comp_run), run_unique_name)

model_stats <- do.call(rbind, 
                       lapply(list.files(paste0("output/update/model_stats/", run_unique_name)
                                         , full.names = TRUE),
                              read.csv))
model_stats <- model_stats[,-1]


print(paste0("CV AUC mean: ", round(mean(model_stats$auc),2)))

## [1] "CV AUC mean: 0.84"

print(paste0("CV RMSE mean: ", round(mean(model_stats$rmse),2)))

## [1] "CV RMSE mean: 0.41"

model_rel_inf <- do.call(rbind, 
                       lapply(list.files(paste0("output/update/model_rel_inf/", run_unique_name)
                                         , full.names = TRUE),
                              read.csv))
model_rel_inf <- model_rel_inf[,-1]

model_rl_long <- model_rel_inf %>% pivot_longer(cols = everything())

means <- model_rl_long %>%
  group_by(name) %>%
  summarise(mean_rl = mean(value)) %>%
  arrange(desc(mean_rl))

model_rl_long$name <- factor(model_rl_long$name, levels = means$name)

ggplot(model_rl_long) +
  geom_jitter(aes(x = name, y = value), size=0.2, width=0.2, alpha=0.5) +
  theme(axis.text.x = element_text(angle = -45,hjust=0)) +
  xlab("Covariate") + ylab("Relative Influence")

ggplot(model_rl_long) +
  geom_jitter(aes(x = name, y = value), size=0.2, width=0.2, alpha=0.5) +
  scale_y_continuous(trans = 'log10')+
  theme(axis.text.x = element_text(angle = -45,hjust=0)) +
  xlab("Covariate") + ylab("Relative Influence (log10 adjusted)")

effect_plot <- readRDS(paste0("output/update/", run_unique_name , "_effect_plot.Rds"))

ggplot(effect_plot) +
  geom_point(aes(x = x, y = value), alpha=0.05, size=0.1) +
  facet_wrap(~var_name)